Theory of Nonlinear Susceptibility and Correlation Length in Glasses and Liquids 
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Within the framework of the effective potential theory of the structural glass transition, we calculate 
for the p-spin model and a hard sphere liquid in the hypernetted chain approximation a static 
nonlinear susceptibility related to a four-point density correlation function, and show that it grows 
and diverges in mean field with exponent 7 = 1/2 as the critical temperature T c is approached from 
below. When T c is approached from above, we calculate for the p-spin model a time dependent 
nonlinear susceptibility and show that there is a characteristic time where this susceptibility has a 
maximum, and that this time grows with decreasing T. We find that this susceptibility diverges as 
T c is approached from above, and has key features in common with the "displacement-displacement 
susceptibility" recently introduced to measure correlated particle motion in simulations of glass- 
forming liquids. 
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Tempted by the possibility of treating the glass tran- 
sition within the framework of conventional critical phe- 
nomena, researchers have long searched for evidence of 
a static correlation length that becomes large as the 
glass transition is approached. However, no clear evi- 
dence for such a length has been reported Q. Recent 
numerical studies of glass-forming liquids have identi- 
fied a dynamical length associated with the range over 
which particle motions are correlated By intro- 

ducing a "displacement-displacement" correlation func- 
tion and generalized "displacement susceptibility" Xu(t), 
Refs. showed in two different model liquids that this 
length grows with decreasing temperature T as the mode 
coupling temperature T c is approached from above, de- 
spite the fact that density and composition correlations 
remain short-ranged. Evidence for a growing length as- 
sociated with correlated particle motion at fixed T in 
simulations below T c has also been reported §. 

In this Letter we calculate within the effective poten- 
tial theory a new susceptibility \i associated with a four- 
point density correlation function jj| . We show that be- 
low T c , X4 diverges with exponent 1/2 as T — > T c 7. We 
show that above T c , \i is time-dependent and resembles 
the susceptibility calculated in [[|||. In particular, we 
show that Xi{t) nas a maximum which diverges in mean 
field as T — > T c + . We argue that (1) the diverging cor- 
relation length implied by the diverging susceptibility is 
associated with incipient ergodicity breaking at T c , and 
(2) this length underlies the growing range of correlated 
particle displacements measured in Refs. |4|||. Finally, 
we test our theoretical predictions above T c using data 
from molecular dynamics simulations of a model glass- 
forming liquid. 

We use two different approaches depending on T. In 



the low T regime (T < T c ), we calculate a static, non- 
linear susceptibility using the effective potential theory 
JlOt , and in the high T regime (T > T c ), we calculate 
a dynamic, nonlinear susceptibility in a dynamical ap- 
proach. The low T calculations are performed both for 
a hard-sphere liquid in the hypernetted chain (HNC) ap- 
proximation [Q and for the spherical p-spin model |Ll] ; 
the dynamical calculations are performed only for the 
spherical p-spin model. This is the simplest model that 
(i) allows both static and dynamic quantities to be calcu- 
lated exactly, and (ii) exhibits several key features com- 
mon to liquids in and close to their glassy regime [^2[dl4| . 
For example, the high-T dynamics of the p-spin model is 
described exactly by the ideal mode coupling equations 
P"5| , |l6[ and displays a sharp transition at T c . 

In the effective potential theory, below T c phase space 
is split into separate ergodic components which remain 
stable at all temperatures (this stability is due to the 
mean field nature of the theory) . A key hypothesis of our 
approach is that in real systems the ergodic components 
are metastable states with a long but finite escape time. 
In this framework we can use a static approach to com- 
pute the properties of these metastable states below T c . 
We can also obtain information on the dynamical prop- 
erties above T c , where mean field theory predicts that 
the time to escape from a metastable state diverges on 
approaching T c . The divergence of a static susceptibility 
inside the metastable state when we approach T c from 
below is directly related to the divergence of the corre- 
sponding dynamical susceptibility when we approach T c 
from above. As is typical for mean field theories, which 
neglect spatial fluctuations, a diverging correlation length 
is deduced from the diverging susceptibility. 

We first describe the essential elements of the effective 
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potential theory (a complete description can be found 
in [jl3||). The theory is formulated using a measure of 
the similarity or "overlap" q between two configurations 
X and Y as an order parameter to detect vitrification. 
Different definitions of q can be used in different systems 
and the main results of the theory do not depend on 
the particular definition adopted. In the case of simple 
liquids with N particles at fixed density |l^Jl^], one can 
define 



1 N 

q(x,Y) = -Vj) 



(i) 



N 



dxdyw(x - y)p x (x)p Y {y), 



where X = {aq, ...,x N }, Y = {yi, ...,y N }, and p z {x) = 

Ylu=i fi( x ~ z i) i s ^ ne microscopic density corresponding 
to the configuration Z = X, Y. Here w(r) is chosen to be 
a smooth, continuous, short-range function close to one 
for r < aro and close to zero otherwise (ro is the radius of 
a particle). The value of a < 1 is arbitrary, and a = 0.3 
is a good compromise for an overlap insensitive to small 
thermal fluctuations uSRi 
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FIG. 1. The effective potential V(q) for the p-spin model, 
for several values of T. At high T the potential is everywhere 
convex, and at low T, V(q) exhibits two minima. In the in- 
set we show the effective potential for a hard sphere fluid 
in the HNC approximation for several values of the density 
(p = 1.0, 1.14, 1.17, 1.19, 1.20). Here the potential is calcu- 
lated around the high- and in the low-q minima, and the lines 
joining the two minima are guides to the eye Jl3|. 

The effective potential V(q), which is a constrained 
free energy, is constructed by choosing a fixed reference 
equilibrium configuration Y at temperature T, and cal- 
culating the free energy of a configuration X that has an 
overlap q with Y: 



V(q) = —log J dXexp(-0H(X))6(q(X,Y)-q). (2) 

Here f3 = 1/fceT and H(X) is the potential energy of X. 
V{q) is self-averaging with respect to the choice of Y. 

The typical mean-field shape of V(q) for a system un- 
dergoing vitrification is shown in Fig. |l| for several values 
of T. The shape of V(q) allows one to distinguish the liq- 
uid from the glassy phase since the presence of a single 
or multiple minima indicates either ergodicity or broken 
crgodicity, respectively. At high T, the system is ergodic 
and V(q) is convex, with a single minimum at a small 
value of the overlap q between any two configurations 
chosen with the Boltzmann weight. Upon lowering T, 
the curvature changes sign, and at T c , V(q) develops a 
secondary minimum at a higher value of q. This signals 
breaking of ergodicity: at T c the configuration space be- 
come disconnected into an exponentially large number 
of "ergodic components" M ~ exp(A^E), each carrying 
vanishing weight in the Boltzmann distribution Jl3|]. A 
fundamental result of the theory, central in the following 
discussion, is that physical quantities calculated in the 
primary minimum represent averages computed with the 
Boltzmann weight, while the same quantities calculated 
in the secondary minimum represent averages computed 
only within a single ergodic component. 

To calculate physical quantities in the effective poten- 
tial theory, we introduce the Legendre transform of V(q): 
r(e) = min g V(q) — eq, where e is a "field" conjugate to 
q, and corresponds to a coupling between configurations. 
For example, the average overlap (q) can be computed 
as (q) = ^| e ^ > where (• • •) represents either of the two 
types of averages. The overlap susceptibility is defined 
as 



X4 



d(q) 



+o =/?iV(V>-<<7> 2 ), 



(3) 



where q = q(X, Y). Inserting Eq. [I] in Eq. || allows us to 
rewrite xa as 



Xi = Nf3 J dx 1 dy 1 dx 2 dy 2 w(x 1 - yi)w(x 2 - y%) 

G i (x 1 ,y 1 ,x 2 ,V2), (4) 

where we have defined the four-point density correlation 
function ||, 

G 4 (xi,yi,x 2 ,y 2 ) = {px{xi)py (vi)px(x 2 )py (2/2)) 

~(px(xi)pY{yi))(px(x 2 )p Y (y2)) ■ (5) 

The two types of averages for Xi are easily calculated. 
We find that when calculated with respect to the Boltz- 
mann average, xa is regular (and small) at all T. How- 
ever, when calculated within the secondary minimum 
(i.e. averaged within a single ergodic component) Xi 
grows for increasing T, and diverges at T c as a power 
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law \i ~ (2c — I 1 ) -7 , as shown in Fig. 2. This demon- 
strates that equilibrium configurations within a single er- 
godic component are highly correlated, while configura- 
tions in different components are not. In both the p-spin 
model and hard sphere model in the HNC approximation, 
the form of V(q) is cubic around the second minimum, 
and thus the value of the exponent 7 is equal to 1/2 
(i.e. the coefficient of the quadratic term resulting when 
V(q) is expanded around the second minimum vanishes 
as (T c — T) 1 / 2 ). This value is universal within mean field 
and coincides with the value of 7 for mean-field spinodal 
transitions. If we assume usual scaling with this value of 
7 then the correlation length exponent v is related to the 
anomalous dimension r\ by v — 1/2(2 — 7]). We empha- 
size that the mean field level at which we are describing 
the system should be the same as that of the ideal mode 
coupling theory. The success of MCT in predicting rela- 
tions between exponents of various dynamical quantities 
leads us to speculate that the mean-field value of 7, 
or a close value, could be observed in real systems, and 
we test this using MD simulations at T > T c later in this 
paper. 
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FIG. 2. The static susceptibility Xi (circles) calculated 
at low temperature (T < T c ), plotted vs. T, together with 
the maximum of the time dependent susceptibility %4(t) (di- 
amonds) at high temperature (T > T c ). The solid line 
through the low-T data is included as a guide to the eye. 
The dashed line through the high-T data indicates a power 
law fit: Xi(T\ = a /( T - T c ) 1/2 + b . In the p-spin model, 
T c = 0.612 [n) and the fit gives a = 3.67, b = -6.28. In real 
systems one can expect a rounding of the divergence at T c . 

We now turn to the high temperature region T > T c 
where the system is ergodic. In this region there is no 
secondary minimum in V(q), and the static susceptibil- 
ity does not exhibit any singular behavior. However, in 
this temperature regime particles of a supercooled liquid 
"oscillate" within cages formed by their neighbors, and 



the system is effectively "frozen" for a characteristic time 
which grows and diverges as T c is approached. This tran- 
sient localization corresponds to highly correlated regions 
of phase space that have finite lifetime and represent the 
high temperature precursors of the low temperature er- 
godic components. 

We study the dynamics of such a system starting in an 
equilibrium initial condition Y = X(0) and evolving in 
time with potential energy H = H[X] — eq(X, Y), and we 
calculate the dynamic susceptibility Xi(t) associated with 
the time dependent overlap q(t) = q(X(t),X(0)) from 
Eqs. 4 and 5, with q = q(t). To calculate Xi(t) above T c 
we use the mode coupling approximation (again the cal- 
culations are performed for the spherical p-spin model, 
where MCT is exact). This gives rise to closed equa- 
tions for the time-dependent overlap correlation and re- 
sponse functions C(t, f) = (q(X(t), X(t' ))) and R(t,t') = 
g^ny which are a slight generalization of the equations 
discussed in Ref. jl0| (the details of the calculation will 
be given elsewhere). We solve for Xi(t) — sc ^'°' > by inte- 
grating these equations numerically as in [ fl0| for different 
values of T. As shown in Fig. 3a, we find that Xi(t) dis- 
plays a maximum as a function of time, which increases 
and shifts to larger t as T — > T+. 
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FIG. 3. (a) The nonlinear susceptibility x^(t) computed 
from the theory in the p-spin model with p = 3, for tem- 
peratures T = 0.7,0.8,0.9,1.0 (T c = .612). The long time 
limit corresponds to the static susceptibility which for this 
model above T c is equal to X4,(°°) = 1/kT. (b) The nonlin- 
ear, time-dependent susceptibility %4(t) calculated for the LJ 
binary mixture above T c , The long time limit for the liquid 
is negligable due to the normalization. Inset: The maximum 
X4.{tt) plotted as a function of T-T c , with T c = 0.435 §0,0. 
The solid line indicates a power law fit to ^4(^4) ~ (T — T c ) , 
and is included in order to compare the simulation data with 
the analytical mean-field prediction. 
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The temperature dependence of the maximum of Xi(f) 
is shown in Fig. 2 (circles); we find that the maximum 
behaves as a negative power of T — T c on approaching T c . 
Although we did not attempt to compute the value of 7 
above T c analytically, in analogy with spinodal points we 
expect that the values of the exponents on the two sides 
of the transition would be the same. Indeed, we find 
that the data are compatible with the mean-field value 
of 1/2 calculated below T c (a power-law fit to the high- 
T data with T c fixed and 7 as a free parameter gives 
7 = 0.52). Thus we find that the dynamic theory of 
the p-spin model predicts a diverging dynamic, nonlin- 
ear susceptibility, (and thus by standard arguments a 
diverging dynamical correlation length), associated with 
a four-point density correlation function. 

In real systems, we expect that the transition at T c will 
be smeared by the existence of dynamical processes that 
restore ergodicity, and which are not taken into account 
in the mean field approach. This has two important ram- 
ifications. The first is that even below T c , Xi{t) should 
display a maximum at finite time. The second is that the 
divergence of Xi as a function of T should be smoothed. 

To test our predictions, we calculate Xi(t) f° r the same 
80:20 Lennard- Jones liquid (containing 8000 particles) 
studied in Ref . 7) . Complete details of the simulation 
may be found in T?]. We evaluate X^if) using the time- 
dependent generalization of Eq. 3, by calculating the fluc- 
tuations in the overlap q measured between two configu- 
rations of the system separated by a time t. In Fig. 3 we 
show Xi{t) as a function of the time interval t between the 
two equilibrium configurations, for seven different values 
of T approaching T c — 0.435 from above. In qualitative 
agreement with our theoretical predictions, we find that 
for the binary Lennard- Jones liquid, x<k{t) nas a maxi- 
mum Xi{tt) a t an intermediate time t\. The amplitude 
of the peak grows and shifts to longer times with decreas- 
ing T. As shown in the inset, the T-dependence of Xii^V) 
is compatible with the mean field prediction. However, 
we caution that a rigorous test of the theory would re- 
quire additional simulations closer to T c and improved 
statistics. 

In summary, we have calculated both within the ef- 
fective potential theory and in a dynamical approach, a 
diverging susceptibility below and above the mode cou- 
pling dynamical critical temperature, respectively. As 
seen clearly from Eq. 4, this susceptibility is related to the 
growing range of a four-point, time-dependent density 
correlation function. Our findings suggest Jh],|21|] an in- 
terpretation of the physics underlying the displacement- 
displacement correlation function calculated in Refs. 
in terms of a four-point density correlation function. The 
correlation function calculated in Refs. ^||| measures the 
extent to which the (scalar) displacements of a pair of 
particles separated by a distance r are spatially corre- 
lated. Specifically, this function is similar to the static, 
two-point pair correlation function g(r), but with each 



particle's contribution to g(r) weighted by its subsequent 
displacement over a time interval [0, f]. In contrast, the 
four-point function G4 studied in the present work mea- 
sures the extent to which "overlapping" particles within 
a time interval [0, t] are correlated. Although these two 
correlation functions are different by definition (and thus, 
e.g. the exponents describing the T-dependence of the 
generalized susceptibilities X4 an d Xu will be different), 
we believe that the length scale associated with G4 un- 
derlies the growing range of correlated particle displace- 
ments measured in Ref. A further critical test of 
our theoretical predictions would consist in the simula- 
tion or experimental measurement of Xiifij or the related 
correlation length below T c . The quantities studied here 
should be measurable experimentally in colloidal liquids 
using confocal video microscopy |22|. 

S.F. would like to acknowledge the kind hospitality of 
the NIST Center for Theoretical and Computational Ma- 
terials Science where part of this work was accomplished. 



[2] 

[3] 
[4] 
[5] 

[6] 

[7] 

[8] 
[9] 



[10] 

[11] 

[12] 
[13] 

[14] 
[15] 



R. L. Leheny, N. Menon, S. R. Nagel, D. L. Price, K. 
Suzuya, and P. Thiyagarajan, J. Chem. Phys. 105, 7783 
(1996); A. van Blaaderen and P. Wiltzius, Science 270, 
1177 (1995). 

Y. Hiwatari and T. Muranaka, J. Non-Cryst. Sol. 235- 
237, 19 (1998); D. Perera and P. Harrowell, ibid, 314.; 

A. Onuki and Y. Yamamoto, ibid, 34. 

B. Doliwa and A. Heuer, Phys. Rev. Lett. 80, 4915 
(1998). 

C. Donati, S. C. Glotzer and P. H. Poole, Phys. Rev. 
Lett. 82, 5064 (1999). 

C. Benneman, C. Donati, J. Baschnagel and S.C. Glotzer, 

Nature, 399, 246 (1999). See also p. 207. 

W.Kob, C.Donati, S.J.Plimpton, P.H. Poole and S.C. 

Glotzer, Phys. Rev. Lett. 79 (1997) 2827. 

C. Donati, S.C. Glotzer, P.H. Poole, W. Kob, and S.J. 

Plimpton, 60, 3170 (1999). 

G. Parisi, J. Phys. Chem, 103(20), 4128 (1999). 
This four-point correlation function was first studied in 
a supercooled liquid in C. Dasgupta, A.V. Indrani, S. 
Ramaswami, and M.K. Phani, Europhys. Lett. 15 (1991) 
307 [Addendum: Europhys. Lett. 15 (1991) 467]. In that 
paper, a negative result was found, possibly because some 
of the data analyzed were in the aging regime, where 
more complex analysis should be done. 
S. Franz and G. Parisi, J.Physique I 5 (1995) 1401; Phys. 
Rev. Lett. 79 (1997) 2486; Physica A 261, 317 (1998). 
A review of the p-s pin model can be found in A. Barrat 
cond-mat/9701031 (unpublished). 

M. Mezard and G. Parisi, J. Phys. A 29 65155 (1996). 
M. Cardenas, S. Franz, G. Parisi, J. Phys. A: Math. Gen. 
31 L163 (1998); J. Chem Phys. 110, 1726 (1 999). 



M. Mezard and G. Parisi cond-mat/9812180 



T.R. Kirkpatrick and P.G. Wolynes, Phys. Rev. A 35, 



4 



3072 (1987); T.R. Kirkpatrick and D. Thirumalai, Phys. 
Rev. B 36, 5388 (1987); T.R. Kirkpatrick and P. G. 
Wolynes, Phys. Rev. B 36, 8552 (1987). 
[16] A. Crisanti, H. Horner and H.J. Sommers, Z. Phys. B 92, 
257 (1993). 

[17] G. Parisi, J. Phys. A: Math. Gen. 30, L765-L770 (1997). 

[18] W. Gotze, in Liquids, Freezing and the Glass Transition, 
J. P. Hansen, D. Levesque, J. Zinn-Justin eds. North Hol- 
land 1990. 

[19] S.C. Glotzer, C. Donati and P.H. Poole, in Computer 
Simulation Studies in Condensed Matter Physics XI 
(Springer- Ver lag, Heidelberg, 1999). 

[20] W. Kob and H.C. Andersen, Transport Theory and Stat. 
Phys. 24(6-8), 1179 (1995). 



[21] S. Franz and G. Parisi, unpublished comment, cond- 



mat 



[22] A.H. Marcus, J. Schofield and SA. Rice, Phys. Rev. E. 
60, 5725 (1999); W. Kegel and A. van Blaaderen, 
preprint; E. Weeks, J.C. Crocker, A.C. Levitt, A. 
Schofield, and D.A. Weitz, preprint. 



5 



